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Premixed turbulent flame calculation 

By S. EL-TAHRY 1 , C. J. RUTLAND 2 , 

J. H. FERZIGER 2 AND M. M. ROGERS 3 

The importance of turbulent premixed flames in a variety of applications has 
led to a substantial amount of effort towards improving our understanding of these 
flames. Although these efforts have increased our understanding, many questions 
remain. For example, there is contradicting evidence concerning the role of the 
turbulent length scales on turbulent flame speed. There are a variety of proposals 
that attempt to correlate the turbulent flame speed with the turbulence velocity 
( e.g ., Bray, 1980 and zur Loye and Bracco, 1987). As in non-reacting turbulent 
flows, the main impediment is the inability to both adequately control experiments 
and to adequately measure the quantities of interest in the experiments. With recent 
advances in supercomputers and the accompanied development of direct numerical 
simulation (DNS) it might now be possible to alleviate some of these difficulties. 

There are a variety of questions that are currently being raised in turbulence 
combustion modeling work for which DNS can provide answers. For example, what 
is the instantaneous structure of a premixed turbulent flame? Flame structures are 
difficult to obtain experimentally because of the fine resolution required. Knowledge 
of the flame structure can aid in closing the viscous dissipation term in the equation 
for the variance of relevant scalars (for a related example see Pope and Anand, 
1985). It is anticipated that in the thin flame mode of combustion ( i.e ., when the 
turbulence strain rate is small compared to the reciprocal of the chemical time scale), 
the flame structure is similar to an undisturbed laminar flame. We need to test this 
assumption under a variety of circumstances such as different chemical systems, 
different Lewis numbers, and different ratios of turbulent velocity to laminar flame 
speeds {u 1 /Si). 

A fundamental question that can be addressed with DNS is what factors control 
the turbulent flame speed (St)- For example, what are the turbulence characteris- 
tics that control St and what is the form of the correlation of these characteristics? 
In DNS we know the values of all quantities at each grid point at every time step. In 
addition, both Eulerian and Lagrangian (i.e., following the flame surface) analyses 
can be conducted and parameters can be switched on and off. Thus with DNS more 
information can be gained than is possible experimentally and many questions can 
be addressed. However, there are limitations to DNS and these are reviewed next. 

The first limitation of all DNS is the range of scales which can be resolved. Small 
grid spacing is required to resolve small scales and large domain size is needed 
to capture the large scales. Computer resources limit this range and hence limit 
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FIGURE 1. Schematic of the flame geometry. 


the maximum Reynolds number of the computations. In the presence of premixed 
flames there is a more stringent condition on the Reynolds number that arises if 
calculations are to be made in the more interesting thin flame regime that occurs 
at high Damkohler numbers (Da). We define Da based on the strain rate of the 
turbulence i.e., 


Da = 


A St 

8u' 




where A and 6 are the Taylor microscale and the flame thickness, respectively. Then, 
it is generally accepted that the thin flame regime is achieved with Da greater than 
unity. To reveal the extra constraint on the Reynolds number in the thin flame 
regime, we can rewrite Da as 


Da = (NiRel^/Nin) 2 , 


( 1 ) 


where Nj is the number of grid lines spanning the computational domain along 
a coordinate direction, N 2 is the number of grid points that will span the flame 
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thickness, Rei is the Reynolds number based on the integral length scale and n is 
the number of integral length scales present in the computational domain. 

In deriving ( 1 ) we used the standard relations between the turbulent scales, hence 
the equation is to be used only for order-of-magnitude purposes, particularly at 
small Reynolds numbers. We notice that N\ is constrained by the maximum job 
size that can be run and N 2 is limited by resolution considerations. We cannot set n 
too low; otherwise the large eddies become too large for the computational domain 
before there has been sufficient time for the flame to develop and propagate. Thus, 
the only avenue available to obtain large Damkohler numbers is to reduce Rei . 

The other relevant constraints on the DNS calculations are the periodic boundary 
conditions and the limitation to incompressible flows. These constraints are par- 
ticular to the computer code that is used in the present study (Rogallo, 1981) and 
in principle can be eliminated. The former constraint may not be too significant. 
The latter constraint, however, precludes studying important aspects, such as the 
influence of the flame on the turbulence and the effect of hydrodynamic flame insta- 
bilities on the flame structure. The present work is an initial effort in studying thin 
premixed flames using DNS. We attempt to address some of the questions raised 
earlier and we do so under the mentioned constraints. 

The problem considered is shown schematically in Figure 1 . It consists of a planar 
flame sheet, located in the y-z plane, initiated at the midpoint of the x-direction. 
The chemical kinetics model used is a hypothetical, single-step kinetics with a single 
reactant A going to a product J3. The equations which govern the thermal and 
chemical states are the conservation equations for species A mass fraction, Y^, and 
the temperature, T. These are: 


D t = 7 y V 2 y„ - bY A exp (~T a /T) 

( 2 ) 

DT 

— = 7 t V 2 T - bHY A exp( — T a /T) , 

(3) 


where D / Dt is the substantial derivative, 7 y and 77 * are the (uniform and constant) 
diffusivities of mass and heat, respectively, b is a pre-exponential factor, T a is the 
activation temperature, and H is the enthalpy of the reaction. In the present 
calculation we have taken the diffusivities of mass and heat to be equal (i.e., the 
Lewis number is unity). For Lewis number not equal to unity, a term involving 
species diffusion should be added to the temperature equation. The value of H was 
selected so as to yield a temperature rise of 15° across the flame. This ensured a 
small density variation due to the reaction so that our assumption of a divergence- 
free velocity field is acceptable. The values of b and T a were selected to give a flame 
thickness spanning about 10 grid lines and an inner flame region spanning about 4 
grid lines. The values of these quantities were established by trial and error using 
a laminar flame code. 

To implement ( 2 ) and (3) in the DNS code (Rogallo, 1981), the source terms 
had to be added to the already existing scalar equations and a chemical time step 
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FIGURE 2. Laminar temperature profiles at 400 step intervals. 
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FIGURE 3. Velocity energy spectra (a) Initial spectra (b) Spectra after 50 steps, 
at the beginning of the reacting part of the calculation. 


stability constraint had to be implemented. The calculations were performed with 
a 128 x 128 x 128 grid. 

To test the accuracy of the flame calculation with the DNS code, a laminar flame 
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propagation calculation was made on the same geometry as depicted in Figure 1 
and the results were compared to the results of a similar calculation made with an 
accurate laminar flame code. The results of the laminar flame calculations made 
with the DNS code are shown in Figure 2. The figure shows the temperature profiles 
at 400 step intervals which begin at t = 0 with a Gaussian distribution. It is seen 
that after 400 time steps the profile has significantly steepened relative to the initial 
conditions and the flame thickness has approximately stabilized (i.e., the initial 
conditions are forgotten). Comparison of the flame speeds obtained from the two 
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FIGURE 5. Temperature contours in x-y planes (a) 2 station = 0 (b) 2 station = 
1/4 (c) 2 station = 1/2 (c) 2 station = 3/4. 


codes were found to be within 2% of each other. This is good agreement considering 
that only about 10 grid points were within the flame and that no de- aliasing was 
attempted for the non-linear source terms. 

For the turbulent combustion calculations, the underlying hydrodynamic field 
was an unforced isotropic turbulence. To achieve a developed turbulence spectrum 
at the start of the chemistry calculations, the latter calculations were initiated after 
developing the hydrodynamics field for 50 time steps. The turbulence spectrum at 
t — 0 and after 50 time steps are shown in Figure 3. The Reynolds number based 
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FIGURE 6. Velocity vectors in x-y planes (a) 2 station = 0 (b) 2 station = 1/4 (c) 
2 station = 1/2 (c) 2 station = 3/4. 


on the Taylor microscale was approximately 5, Da was approximately equal to 1.5 
and u' /Si was of the order 0.53. This meant that we expect a thin flame with low 
turbulence (i.e., the wrinkled flame regime). 

We start reviewing the results by examining the flame structure. Figure 4 shows 
the constant temperature surfaces for T = 302° and T = 314° (the unburned 
temperature is 300° and the burned temperature is 315°). The surfaces are seen 
to have a wavy topology consistent with what is expected in the wrinkled flame 
regime. Figure 5 shows contours of constant temperatures in x-y planes taken at 
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FIGURE 7. PDF of the mean reactant mass fraction. 


the 0, 1/4, 1/2, and 3/4 stations in the ^-direction. Here the wrinkled structure is 
more evident and the wrinkling seems to be occurring at several scales. The flame 
thickness is found to vary with location. This variation is caused by the inclination 
of the local flame front relative to the section plane. Thus the flame appears thicker 
than it is along a normal to the flame front. The flame thickness in the thin regions 
were found to be comparable to the thickness of the laminar flame calculated earlier. 
At thirteen locations tested the flame structure was also found to be similar to the 
laminar flame structure. This result agrees with the first-order solution derived 
by Clavin and Williams (1979). However, their second-order solution modifies the 
laminar profile. 

Figure 6 shows the velocity vector field in x-y planes at the same z- locations as 
Figure 5. In these figures the flame wrinkling appears for the most part to follow 
the velocity field. We can also see that on the scale of the flame thickness the strain 
rate is small. This qualitatively verifies that we are in the thin flame regime. 

The probability density function (PDF) for the reactant mass fraction is shown 
in Figure 7 at various x-locations. These PDFs were generated by collecting data 
from the homogeneous planes (y-z planes) at the specified x locations, dividing it 
into 20 bins and averaging the data in each bin. The distributions are seen to be 
qualitatively similar to the PDF of the laminar flamelet model proposed by Bray 
and Moss (1977). According to this model the PDF consists of two delta functions 
located at Ya = 0 and Ya = 1 which are joined by a region of low probability 
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FIGURE 8. Various scales and velocities as a function of time. □ Taylor microscale; 
o u'; ^ Reynolds number/10; + Turbulent flame speed; Laminar flame speed. 


that has its PDF proportional to the inverse of the gradient of Y A . We notice from 
Figure 7 that, because the gradients on the cool side of the flame are small, the 
probability density function tends to be higher at the large values of Y A . There 
is, however, a peculiar behavior of the PDFs at some of the x-locations where 
we observe a maximum of the PDF at values of Y A of around 0.9. This probably 
indicates that the wrinkling of the flame is not much larger than the flame thickness. 
This is supported by the contour plots in Figure 5. 

Next we will examine the dependence of the turbulent flame speed on u'. Figure 8 
depicts the variation of A,u', Re\ and the turbulent flame speed St with time. The 
laminar flame speed is also included for comparison. The hydrodynamic quantities 
display the expected behavior i.e., u' and Re\ decrease with time, while A increases. 
Hence, Da will increase with time which means that we remain in the wrinkled 
flame regime for the duration of the calculation. The variation of the turbulent 
flame speed with u' is shown in Figure 9. We find St increasing with u' but the 
increase is much slower than suggested by the commonly used formula 


St = Si + c'u' , 


(4) 
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FIGURE 9. Flame speed vs. u'. □ Turbulent flame speed; Williams’ equation; 

Laminar flame speed. 


where c' is a constant. Tabaczyski et al. (1974) use c' = 1 , Witze and Mendes- 
Lopez (1986) suggest a value of 2.2. For other suggestions see zur Loye and Bracco 
(1987). Our results are closer to the theoretically obtained formula given by Clavin 
and Williams (see Williams 1985) for weak turbulence 


St = Si 


^(1 + y/\ + 8c(u'/ S l) 2 ) 


( 5 ) 


where c is a constant of order unity. This equation is also plotted in Figure 9 using 
c = 1.2 to match the first data point. Within the range of our data the magnitude 
of St/Sl given by equation (5) is fairly close to our calculated values. However, 
the trend is not, indicating that St probably depends on more than just v' . It is 
interesting to note that in equations (4) and (5) and most other suggestions in the 
literature, there is no dependence of St/ S i on heat release or Reynolds number 
(there are some exceptions e.g ., Abdel-Gayad and Bradley (1976) give an expression 
that involves the Reynolds number). Hence the flame speeds calculated with the 
DNS code should also be representative of flames with large density variation and 
large Reynolds numbers. However, as implied earlier, the differences in the various 
proposed equations for the flame speed may be due to not including all of the 
relevant parameters on which the flame speed depends. We should try to investigate 
these systematically. 

Finally, the development of the mean temperature profiles is shown in Figure 10. 
The profiles suggest that the turbulent flame thickness is about twice the laminar 
flame thickness. We also observe that the flame thickness increases with flame 
travel. This may be due to the increase of the turbulent length scale with time. 
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FIGURE 10. Mean turbulent temperature profiles at 200 step intervals. 


Our plans to continue the work on this project include more analysis of the case we 
have presented here. In particular, we will examine scales and structures relating to 
the turbulent diffusion of the reactant mass fraction. We will also calculate the area 
of the wrinkled flame front and compare this to the turbulent flame speed. Further 
work in this area will include turbulent calculations with a mean strain-rate so that 
we can examine the effects of anisotropy on flame propagation. 
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